
#ifndef AFTERMATH_RANDOM_ZIGGURAT_CORE_HPP_INCLUDED
#define AFTERMATH_RANDOM_ZIGGURAT_CORE_HPP_INCLUDED

#include "../template_math.hpp"

#include <cmath>
#include <cstdint>
#include <functional>
#include <type_traits>

namespace aftermath
{
    namespace random
    {
        namespace detail
        {
            template <std::size_t t_n_boxes, typename t_uniform_type>
            struct ziggurat_choose_box
            {
                typedef t_uniform_type uniform_type;

                static std::uint_fast32_t choose_box(uniform_type uniform_random);
            };

            template <std::size_t t_n_boxes>
            struct ziggurat_choose_box<t_n_boxes, std::uint_fast32_t>
            {
                typedef std::uint_fast32_t uniform_type;

                static std::uint_fast32_t choose_box(uniform_type uniform_random)
                {
                    static_assert(is_power_of_two<t_n_boxes>::value, "<t_n_boxes> has to be a power of two.");
                    return uniform_random & (t_n_boxes - 1);
                }
            };

            //template <std::size_t t_n_boxes>
            //struct ziggurat_choose_box<t_n_boxes, double>
            //{
            //    typedef double uniform_type;

            //    static std::uint_fast32_t choose_box(uniform_type uniform_random)
            //    {
            //        static_assert(is_power_of_two<t_n_boxes>::value, "<t_n_boxes> has to be a power of two.");
            //        static_assert(sizeof(uniform_type) >= sizeof(std::uint_fast32_t), "size of <t_uniform_type> has to be greater than <std::uint_fast32_t>.");

            //        // Take 4 least significant bytes of the mantissa.
            //        std::uint_fast32_t* address = reinterpret_cast<std::uint_fast32_t&>((&uniform_random) + (sizeof(uniform_type) - sizeof(std::uint_fast32_t)));
            //        std::uint_fast32_t discrete = *address;

            //        return discrete & (t_n_boxes - 1);
            //        // Alterntaively, can do multiplication and floors, but the above approach should be much faster.
            //    }
            //};
        }

        template <std::size_t t_n_boxes, typename t_ziggurat_impl_type, typename t_uniform_type, typename t_bounds_type, t_bounds_type t_diameter, typename t_result_type>
        struct ziggurat_core
        {
            static const std::size_t n_boxes = t_n_boxes;
            static const t_bounds_type diameter = t_diameter;

            typedef t_ziggurat_impl_type ziggurat_impl_type;
            typedef t_uniform_type       uniform_type;
            typedef t_result_type        result_type;
            typedef t_bounds_type        bounds_type;
            typedef ziggurat_core<n_boxes, ziggurat_impl_type, uniform_type, bounds_type, diameter, result_type> type;

            template <typename t_engine_type>
            result_type sample(t_engine_type& uniform_generator)
            {
                static_assert(std::is_same<typename t_engine_type::result_type, uniform_type>::value, "type mismatch");
                static_assert(t_engine_type::max() - t_engine_type::min() == type::diameter, "<t_engine_type>::max() - <t_engine_type>::min() has to be equal to <diameter>.");

                auto that = static_cast<t_ziggurat_impl_type*>(this);
                for (;;) // while (true)
                {
                    uniform_type u1 = uniform_generator() - t_engine_type::min(); // Used for box index.
                    uniform_type u2 = uniform_generator() - t_engine_type::min(); // Used to generate the horizontal component.

                    // ~~ Pick a box ~~
                    std::uint_fast32_t box_index = detail::ziggurat_choose_box<n_boxes, uniform_type>::choose_box(u1);

                    // ~~ Optimized rejection ~~
                    bool is_interior = false;
                    result_type z = that->sample_from_box_horizontal(box_index, u2, std::ref(is_interior));
                    if (is_interior) return z;

                    // ~~ Tails ~~
                    if (that->is_tail_box(box_index)) return that->sample_tail(box_index, uniform_generator);

                    // ~~ Standard rejection ~~
                    uniform_type u3 = uniform_generator() - t_engine_type::min(); // Used to generate the vertical component.
                    if (that->is_inside_box_vertical(box_index, z, u3)) return z;
                }
            }
        };
    }
}

#endif // AFTERMATH_RANDOM_ZIGGURAT_CORE_HPP_INCLUDED
